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(A) INTRODUCTION 

As the result of economic and population growth, the demand 
for electrical power has been doubling every decade since 
World War II in the U.S. The rapid increase both in elec- 
tricity production and in the unit size of thermal power 
plants has resulted in use of large quantities of water for 
cooling purposes. The most common system for heated efflu- 
ents dsposal is the open cycle system where the condenser 
cooling water is discharged back into the natural water- 
way. Analysis of thermal discharge with a view to decreas- 
ing environmental effects and increasing the cooling 
efficiency has become an important part of power plant 
siting evaluations. 

In the past years, numerical models have been developed by 
many researchers for predicting the velocity and tempera- 
ture distribution in thermal plumes, such as Wada (1969) , 

Koh and Fan (1971), Stolzenbach and Harleman ( 1972) ,Waldrop 
and Farmer (1974), Paul and Lick (1973) etc. A review of 
existing models is presented by Dunn and Policastro (1975) . 
Most of the existing models cannot adequately account for 
time-dependent, three-dimensional flow fields, realistic 
bottom topography effects, and surface wave phenomenon. 

Since a buoyant plume is three-dimensional in character 
and is significantly effected by bottom topography, ambient 
currents and local meteorology, the need for "complete 
models" is acute. 

A time-dependent, three-dimensional, free-surface , niamerical 
model with both horizontal and vertical stretching has been 
developed for the thermal discharge study, and is presented 
in this paper. The effects of variable bottom topography, 
free surface elevations, surface heat transfer, currents, 
and meteorological conditions are considered in this model. 
The model can simulate the physical conditions needed to 
provide sufficient three-dimensional information for the 
studies of water quality to meet the appropriate standards. 
The model is not only a flexible and economical tool for 
the heated water disposal system design or improvement 
consideration, but also for monitoring thermal effects 
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in the receiving domain. 

Numerical integration of the governing equations (conserva- 
tion of mass, conservation of momentum, and energy) is 
performed simultaneously by using finite difference approxi- 
mations for the differential eq’iation. Both surface and 
submerged types of discharge systems in a bay and coastal 
site have been modelled. The results from the computer 
model were compared satisfactorily both with airborne 
thermal scanner remote sensing data and in-situe ground 
truth measurements at Florida Power & Light Company's 
Cutler Ridge Plant and Hutchinson Island Nuclear Power 
Plant. Only the results for Hutchinson Island site will 
be presented here. The application of this model to a 
specific site, involves specification of initial conditions, 
boundary conditions, and information geometry of the discharge. 

(B) FORMULATION AND BASIC EQUATIONS 

Numerical modeling of heated water discharges involves the 
application of fundamental hydrodynamic and thermodynamic 
equations. Several assumptions have been made for this 
study, namely, 

(1) The fluid flow is incompressible, 

(2) The fluid flow is turbulent and that the molecular 
transport can be neglected in comparison with turbu- 
lent transport, approxi -kited by eddy transport co- 
efficients, 

(3) Pressures are hydrostatic throughout the whole domain, 

(4) The solid wall and bottom are adiabatic, and 

(5) The density throughout tlie flow field is a function of 
temperature only, 

A vertical coordinate transformation techniqiie similar to 
Phillips (1957) is used to incorporate both free surface and 
variable bottom topography and the same number of grid points 
in vertical direction can be used in the model. The new 
sigma vertical coordinate system is obtained by letting 

a = g (x,v, 2 ,t) ^ z + T\ (x,v,t) I 

H (x,y,t) h (x,y) + T) (x,y,t) 

where Z = Z (x,y,z,t) is the position of the fluid element 
relative to the free surface, H = H (x,y,t) is the depth con- 
tour relative to the free surface, z represents the vertical 
position relative to the mean water level, T1 is the free 
surface elevation measured positively upward from mean water 
level, h represents the depth relative to the mean water 
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level. The value of a ranges from zero at free surface to 
unity at the bottom of the basin. This mathematical tech- 
nique has been used for a variety of numerical models of 
geophysical fluid flow, e.g., Smagorinsky (1965), Arakawa 
(1972), for the atmosphere study and Freeman (1972), 

Sengupta and Lick (1974) for the study of lakes. Fig. 1 
and Fig. 2 show the two coordinate systems before and 
after vertical coordinate transformation respectively. 

It is desireable to obtain a more detailed description of 
the flow near the discharging area and meanwhile a large 
horizontal domain can be covered. If a constant grid size 
were used, then a large number of grid points and excessive 
computation time are needed for the model. In order to 
solve this problem, a horizontal stretching is used in both 
lateral and transverse directions to create a small grid 
size near the discharge area and larger grid size away 
from the discharge points. An arc-tangent equation was 
used by Waldrop and Farmer (1974) as the horizontal stretch- 
ing. A hyperbolic sine stretching has been used in this 

study for both lateral and transverse directions, as shown 
below: 


X = a 


+ sinh 


C C„ (X - d) ] 


= b + sinh 


C (Y - e) ] 


where x, y are the real coordinate; X,Y are the stretched 
coordinate; a and b are the distance at which the minimum 
step size is desired; C., C_, C_, C., d,e, are the constants 
to be determined by the'^imposed'^conditions . By using this 
horizontal stretching, even increments of ^ and AY can be 
used in the model for the variable physical grid sizes. 


Applying the differential transformation relationships, the 
horizontal and vertical stretched governing equations in the 
XYa coordinate are; 

Continuity Equation 

X' Y' AiHvJ. ^ H an ^ Q 

3X qY ^ 

Momentum Equation 
u - momentum; 


?^(Hu) X' (Huu) 
ht ^ 9X 


Y' 


+ H 
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hC-^(||) +gx* + fv] 

P 5 X 


+ Kjj [ ( X* )^ || ^ + H ( X' ^ + HX' 


^ -| 


KjjC (Y*)^-^ ^ + H(Y')^ ^ + HY" ^ ] 


by 


by 


^ C nt? < » J'v ) 3 


( 2 ) 


V - momentum; 


a(Hv) ^ ,fl( Huy) , >(Hw) ^ j- ;^(vQ) ^ 

Bt aX ^Y ^CT 


H [ - f- ( || ) + gy ( a |S - ) - f u : 


+ Kjj [ ( X' )^ || ^ + H ( X* )^ -^ + HX" -^ ] 


bx 


bx 


+ Kjj[ (Y' (Y‘ + HY" ^ ] 


+ 1 ( p K ] 


(3) 


Energy Eq uation 


> (HT) + X' ;^(HuT) + Y' ^iHuT) + H ^(OT) 
Bt 9X gY 9ff 


= B„ C ( X' + H ( X' ^ + 


H 


bx 


HX' 


BX‘ 


bx 


+ B„[ (Y* )^-^^ + H (Y’ )^-^ + 


H 


by 9Y 


HY' 


by 


^ ] 

BY 


+ 1 


( Q B .SS. ) *1 

pH ^ ga ' P ®v aa ^ 


ri7 [ 


(4) 
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Eauation of State 


p = F (T) = 1.000428 - 0.000019T - 0o000036T^ (5) 

The relationship between real vertical velocity and vertical 
velocity in cr coordinate is derived as follows: 

w = nn + ^ 


By integrating the continuity equation with respect to a 
from surface to bottom we can get 


9t 


r X* 1 dff 

Jo ^ ax aY ^ 


( w^- X- 


2ih _ 



^ \ 


0 .( 7 ) 


Substitute the abo’»'e equation into continuity equation and 
integrate from surface to cr = a, and the result can be ob- 
tained as 


n 



C X- + Y- 2^ ]a<T + f j-^ 


[ X' 


AiHal 

9X 


+ 


Y' 


a(Hv) 

by 


]dff 


2. ( >^1. 
H 


X' 


ak _ V, , ^ ) 
bx ° aY 


( 8 ) 


where subscript b denotes at the bottom of basin and 


X- = as ; y. = 4Y 


dx 


dy 


... _ ^ 


X" = 


... _ d__Y 


Y" = 


dx 


dy 


The hydrostatic relationship is used as the diagnostic 
equation for pressure 

P (a) = P (0) + gH P (a) dcT (9) 

where P (0) is the pressure at the free surface. 
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(C) INITIAL AND BOUNDARY CONDITIONS 

In order to solve the set of governing equations, the initial 
and boundary conditions should be specified throughout the 
whole domain. This information may be obtained from remote 
sensing data, in-situ measurement and other sources, Bomdary 
conditions are specified at the air-water interface, the 
bottom of basin, the lateral solid walls, discharge points, 
and open boundaries. At the air-water interface the wind 
stress and the heat transfer coefficient are specified. The 
condition of no slip at lateral solid wall and bottom of 
basin is used for the momentum equations. 

At solid lateral walls and at the bottom of basin, adiabatic 
conditions are assumed. At points of thermal discharge, both 
velocity and temperature are specified. The domain is bound 
on three sides by open boxindaries where conditions are most 
difficult to specify. At these boundaries, the outflow 
velocity gradient normal to boiindary is considered zero and 
the inflow velocities are the known values. For the tempera- 
ture, the second derivative is equated to zero for outflow, 
implying no diffusive heat transfer and convective heat trans- 
fer is dominated. For the inflow, the temperature is the same 
as ambient water temperature, 

(D) METHOD OF SOLUTION 

The set of governing equations are unsteady, non-linear and 
second order partial differential equations. With the wide- 
spread use of the high-speed digital computers, finite differ- 
ence approximation methods for the partial differentiation are 
used for solving the numerical solution. 

In the three-dimensional free surface model, central differenc 
ing is used for the space derivatives in the interior of the 
basin. Single-sided differencing is used for the space derive 
tives at the boundary. Forward differencing in time for the 
first time step and then the central differencing in time is 
used for predicting the new values of all the dependent 
variables. In order to insure the numerical stability, the 
convective terms, pressure term and Coriolis force term are 
calculated by using the values at present time, the viscous 
terms are calculated by using the values at previous time step 
(t-At) , All the dependent variables are computed at integral 
grid points of the basin, 

A two-dimensional matrix (MAR) was created in the computation 
process for locating the position in the basin to insure the 
appropriate numerical scheme and boxindary conditions are used 
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in the computation. Different values of MAR were assigned to 
the different locations xn the basin, such as MAR = 11 rep- 

rents the points in the interior, and MAR = 0 represents points 
outside the domain. 

The finite difference form of governing equations is solved 
simultaneously for each time step. Before the computation 
process, the initial conditions for all the dependent 
variables should be specified throughout the domain. The 
surface heights are computed from the vertically integrated 
continuity equation. The horizontal velocities are computed 
directly from the momentum equation. The vertical velocities 
are obtained from new horizontal velocities and surface 
heights into equation (6) , The temperatures are calculated 
from the energy equation. The pressures are obtained from 
the hydrostatic relationship. This procedure is performed 
at each time step until the desired length of time has 
been achieved. 


(E) APPLICATION 

The three-dimensional free surface model has been applied to 
FPL's Hutchinson Island Nuclear Power Plant which is located 
about midway between the cities of Fort Pierce and Stuart 
on the eastern ocean coast of Florida, The objective of this 
study is to use this model to predict the velocity and 
temperature distribution caused by a submerged thermal dis- 
charge into the ocean vinder various environmental conditions. 
The 12 feet diameter submerged discharge pipeline is buried 
in the ocean bed and terminates at a point about 1200 feet 
offshore at a depth of 18 feet from mean water level. At 
its termination, a two port Y-type discharge is added with 
each arm being 7.5 feet in diameter. A short sloping con- 
crete pan is located at the outlet to prevent scouring of the 
ocean floor. The submerged jet discharge exits horizontally 
with a relatively high velocity and rises to the surface 
owing to buoyancy effects and then spreads out by the action 
of turbulent mixing and cimbient currents, 

A numerical grid system of 20 x 20 in the horizontal plane 
and five grids in the vertical plane were chosen. The domain 
covered is 2000 meters parallel to the shore line and 2380 
meters perpendicular to the shore (along the discharge pipe 
axis) , Fig. 4 shows a vertical section along the discharge 
pipe axis, the vertical and horizontal stretching and its 
influence on the nvunerical grid system is apparent. 

On June 2, 1976, the heated water had an exit velocity of 280 
cm/sec (9.1 ft/sec) at each end of the Y-type discharging 
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pipe. The anibient water temperature was about 25.5 C, dis- 
charge water temperature was 35°C, and air temperature was 
29°Co The wind speed was 10 mph from the southeast. The 
ocean current was about 25 cm/sec predominantly northerly. 

The discharge conditions and the environmental conditions 
were incorporated into the model. It was found that the 
plume rapidly rises to the surface and spreads to the north 
owing to the current. 

Fig. 5 shows the surface velocities predicted for environ- 
mental conditions of June 2, 1976. 'The imposed northerly 
current prevails away from the discharge and near the dis- 
charge points there is a source -like flow pattern because the 
plume rapidly rises to the surface cind spreads. Fig. 6 shows 
the horizontal velocity on the plane 6.5 meters below the mean 
water level. The Y-type discharge pipe d4scharged horizontally 
with a relatively high velocity as can be seen at the dis- 
charge points, and then turned with the current. Fig. 7 shows 
the velocity distributions of the vertical section along the 
discharge pipe axis. The vortex just west of the discharge 
caused by the entrainment is clearly evident. Fig. 8 shows 
the velocity distributions of the vertical section I = 10 that 
is parallel to the shore line. The plume rises to the sur- 
face owing to the buoyancy and then is carried along with 
the current. 

Fig. 9 shows the isotherms of a vertical section 1=9 that 
is perpendicular to the axis of the discharge pipe. The 
plume rxses to the surface and then spreads out and is 
carried along with the current as is evident. Fig. 10 shows 
the surface isotherms from I.R. data in the morning at 
Hutchinson Island Power Plant on June 2, 1976. The submerged 
plume reises to the surface and spreads with the northerly 
current. The morning I.R. data base was used as the initial 
conditions for the model and the afternoon results were pre- 
dicted by the model. The plume shifted eastwards during this 
period and the area of 26.5°C isotherm was reduced. The 
width of the thermal plume increased somewhat. Fig. 11 shows 
the comparison of surface isotherms obtained by the model 
prediction and I.R. measurement data. In general, the model 
prediction is observed to be in relatively good agreement 
with I.R. measurements. 

(F) SUMMARY 

The present study predicts the three-dimensional flow and 
temperature field in the receiving body of ocean water owing 
to the submerged discharge near the ocean floor. The effects 
of currents, wind, surface cooling, and bottom topography 
are combined in this complete field numerical model. The 
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vertical and horizontal stretching provide more effxcient use 
of this model and ease of boundary condition specification. 
The results from the model obtained indicate proper numerical 
behavior of the model. It can be concluded that the physical 
behavior of the pl\ame is reproduced by the model with rela- 
tively good agreement. Further calbiration efforts are still 
continuing. Some field experimental data will be obtained 
for using the final verification of the model. 
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Fig. (5^ Surface velocity distribution with current, 
wind and bottom topography at Hutchinson 
Island Site (Free Surface Model) 
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Fig. (6^ Horizontal velocity distribution at the plant 

6.5 meter below the free surface with current, 
wind and bottom topography at Hutchinson 
Island Site (Free Surface Model) 
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Fig. (7 } Velocity distribution at section J-ll for Hutchinson Island Site 
with vertical scale exaggerated 5.25 times 
(Free Surface Model) 
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Fig^(8 J '* Velocity distribution at section I— 10 for Hutchinson Island Site 
with vertical scale exaggerated 5.25 times (Free Surface Model) 
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Vertical isotherms at section 1=9 
for the Hutchinson Island Site 
(Free Surface Model) 
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Hutchinson Island Power Plant, 2 June 
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